function [alphav,betav]=msdtoab(MUvec,SDvec,densvec,lowerb,upperb)
% ================================================================= 
%   function [alphav,betav]=msdtoab(MUvec,SDvec,densvec)
%
%   Find alpha, beta to match the mean and standard deviation of
%   of the densities specified by MUvec SDvec and desnvec
%   Output
%   ======
%   alphav & betav which contain the alpha & beta for each density 
%   For Normal alpha=Mean beta = SDV 
%
%   For W (IW)    following the notation in Lubrano et. al 
%            alpha = d (deg. freedom)
%            beta= = v ( scale ) 
%   such that if sigma^2 ~ IW(v,d) then 1/sigma^2 ~ Gamma( d/2,2/v) 
%
%   For the uniform MU=a < SD=b 
%
%   Created 6/18/03 
%   It is the shorter version of alphabeta.m to use within a program 
% 
%   5/5/04    Added the IG1 distribution using function in 
%             Dynare  
%
%   22/3/06  Bounds of the uniform obtained from the MEAN and SD 
%            Converting into A and B. 
%   1/13/08 Moved inside here the Uniform instead of the adjustment in
%   LOADPRIOR.m 
%
%   1/13/08  Added the T distribution, which a BETA with exact truncation 
% ***********************************************************************
% Storage Matrices
MUvec = MUvec(:); 
SDvec = SDvec(:); 
np    =size(MUvec,1);
alphav=zeros(np,1);
betav =zeros(np,1);
densvec=char(densvec); 
for ii=1:np;
    dens=densvec(ii);
    MU=MUvec(ii);
    SD=SDvec(ii);
    if strmatch(dens,'N');
        alp=MU;
        bt=SD;
    elseif strmatch(dens,'G');
        bt=(SD^2)/MU;
        alp=MU/bt;
    elseif strmatch(dens,'B');
        bt=MU*( (1-MU)^2 )/(SD^2) - ( 1 - MU );
        alp=MU*bt/(1-MU);
    elseif strmatch(dens,'W');
        if SD > 100
            alp = 4;
            bt = MU*(alp - 2);
            fpraid('Treating the variance of the IW as inf');
        else
            alp=2*(MU^2)/(SD^2) + 4;
            bt=MU*(alp-2);
        end
    elseif strmatch(dens,'I');
        if SD > 100 ;
            SD=inf;
            disp('Treating the Variance of IG1 as Infinity');
        end
        [bt,alp]=inverse_gamma_specification( MU , SD ,1 );
    elseif strmatch(dens,'U');
        %[alp,bt]= u_msd2ab(MU,SD,[]);
        alp=lowerb(ii); 
        bt =upperb(ii); 
        if lowerb(ii) > upperb(ii); error('Lower bound greater than Upper bound'); end 
    elseif strmatch(dens,'T'); 
        %disp('Not ready....'); 
        if lowerb(ii) > upperb(ii); error('Lower bound greater than Upper bound'); end
        [alp,bt]=msdbetab(MU,SD,lowerb(ii),upperb(ii)); 
        disp('Using Truncated Beta'); 
        dispaj('MU=',MU,' SD=',SD,' LowerB',lowerb(ii),' UpperB',upperb(ii)); 
        dispaj('Implied ALPHA=',alp,' BETA=',bt);
    elseif strmatch(dens,'E'); 
        
        % Solve for the "k" parameter in weibull
        weibFunc = @(k) gamma(1+1./k).^2./gamma(1+2./k)-MU^2/(SD^2+MU^2);
        grid=(1E-3:.1:100);
        
        % Find an approximate solution.  Note that weibFunc is strictly
        % montone increasing so this rough minumum will be unique
        [val loc] = min(abs(feval(weibFunc,grid)));
        
        % Find a more exact solution
        [bt fval exit_flag]=fzero(weibFunc,grid(loc));
        
        if exit_flag ~= 1
            error('Failed to solve for weibull parameters');
        end
        
        % Solve for "\lambda" parameter
        alp = MU/gamma(1+1/bt);
        
        % Check to see if solved parameters are reasonable
        if (alp < 0) || (bt < 0)
            error('Weibull parameters not in parameter space')
        end
        
    else
        error('Distribution not matched')
    end
    alphav(ii)=alp;
    betav(ii)=bt;
end;